%ֻ������ֵ��Ϊ������������
%clc;
%clear all;   
%iII=1;
function kmeans(num)

iI=num;
t=5;
file='test';
file1=strcat(file,num2str(iI));
p1=strcat(file1,'.tif');

file2=strcat(file,num2str(iI+1));
p2=strcat(file2,'.tif');
t_kmeans1=cputime;
%�任��Сsize350*366
I0=imread(p1);
[Hh,Ww]=size(I0);
% figure,imshow(I0);
%�Աȶ���ǿ
I=ctrenhance(p1);
%I=I0;
%I=double(I);
x21=mat2gray(I);%ת��Ϊ�Ҷ�ͼ�� 0-1
[K11,K21]=size(x21);

%%%%%%%%%%%%%%%%��ֵ�˲��õ�y1%%%%%%%%%%%%%%%%%%
m=3;
I11=medfilt2(x21,[m,m]); 
% figure,imshow(I11)
%��I��������λ��para1����һЩ��������morphspotgrid1��ע��
%para2Ϊ����ˮƽ������λ�õľ���
%para3Ϊ���ش�ֱ������λ�õľ���
%hnum=para1(5);%ˮƽ������
%vnum=para1(6);%��ֱ������
%[para1,para2,para3]=morphspotgrid1(I);

[para1,para2,para3,para4,para5]=morphspotgrid(I11);
para2=round(para2);
if para2(1)==0
    para2(1)=1;
end
para3=round(para3);
if para3(1)==0
    para3(1)=1;
end

hnum=size(para2,2);%ˮƽ������
vnum=size(para3,2);%��ֱ������
%para4
%para5
I=uint8(I11.*256);
[h,w]=size(I);
h1=size(para4,2);
v1=size(para5,2);

% hnum=hnum-3;
% h1=h1-5;
% para4=para4(1:h1);
%para4(h1)=h;
% hnum=hnum-1;
% h1=h1-2;
% para4(1:11)=para41(1:11);
% para4(12:h1)=para41(14:h1+2);

t_kmeans01=cputime;
% t=3;


for j=1:h
    if j<para4(2)-t
        I(j,:)=0;
    end
    for i=3:2:h1-1
       if j>para4(i)+t&&j<para4(i+1)-t
         I(j,:)=0;
       end
    end
    if j>para4(h1)+t
        I(j,:)=0;
    end
end


for j=1:w
    if j<para5(2)-t
        I(:,j)=0;
    end
    for i=1:2:v1-1
      if j>para5(i)+t&&j<para5(i+1)-t
        I(:,j)=0;
      end
    end
    if j>para5(v1)+t
        I(:,j)=0;
    end
end

III=zeros(h,w);
% III2=zeros(h,w);
for i1=1:hnum-1
    for j1=1:vnum-1
        for i2=1:para4(2*i1+1)-para4(2*i1)+2*t-1
            for j2=1:para5(2*j1+1)-para5(2*j1)+2*t-1
                I1(i2,j2)=I(i2+para4(2*i1)-t,j2+para5(2*j1)-t);
%         for i2=1:para2(i1+1)-para2(i1)+1  
%             for j2=1:para3(j1+1)-para3(j1)+1
%                 I1(i2,j2)=I(i2+para2(i1)-1,j2+para3(j1)-1);        
            end
        end
        [H,W]=size(I1);
        
ii=0;
num=H*W;
II=zeros(num,1);
%����������ֵ����������II
for i=1:H
    for j=1:W
        ii=ii+1;
        II(ii,1)=I1(i,j);
    end
end
clear I1


%ʹ��matlab�Դ���k-means����
%[idx,C,sumD,D]=kmeans(II,2,'Replicates',10,'Display','final','MaxIter',1000,'EmptyAction', 'singleton','OnlinePhase','off');
% [idx,C,sumD,D]=kmeans(II,2,'start','uniform');
idx=kmeans2(II,2,20);
idx=2-idx;

%����kmeans�õ����б����idx�б�ǰ���������
iii=0;
for i3=para4(2*i1)-t+1:para4(2*i1+1)+t-1
    for j3=para5(2*j1)-t+1:para5(2*j1+1)+t-1
% for i3=para2(i1):para2(i1+1)
%    for j3=para3(j1):para3(j1+1)
        iii=iii+1;
        %if sum(idx)<num/2-4%�����ٵĵ�Ϊ������������ĵ�Ϊǰ��
%         if idx(1,1)==1
            if idx(iii)==1
              III(i3,j3)=0;
            else III(i3,j3)=1;
            end
%         else
%             if idx(iii)==1
%               III(i3,j3)=1;
%             else III(i3,j3)=0;
%             end
%         end
    end
end
    end
end

% figure,imshow(III)
% imwrite(III,'kmeans_cy1.tif','tif');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�󲹳�����ԭͼ�Ͻ��ָ��Ե����
%kmeans����󻮷ֱ�Ե
kmeans_edge3=edge(III,'canny');
PP3=imread(p1);
[h1,w1]=size(PP3);
figure,imshow(PP3,[]);
for i=1:h1
  for j=1:w1
         if kmeans_edge3(i,j)==1
            hold on
         plot(j,i,'r');
         end
   end
end
% title('kmeans edge of p1 ');

% figure,imshow(PP3,[]);
% hold on
% plot(1:w,para4(2)-4,'r');
% hold on
% plot(1:w,para4(3)+4,'r');
% hold on
% plot(para5(2)-4,1:h,'r');
% hold on
% plot(para5(3)+4,1:h,'r');
time_kmean1=cputime-t_kmeans1
time_kmean01=cputime-t_kmeans01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�������յõ��ķָ�ͼ������������ݵĸ�������
%������������ڵĸ�������
PP1=imread(p1);
PP1=double(PP1);
PP1=PP1./256;        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ǰ���������صĻҶ�ֵ��CH1_RAW�ͱ����������صĻҶ�ֵ��CH1_BKD_RAW
%��ǰ���������ҶȾ�ֵ����CH1_MEAN-CH1_BKD_MEAN=CH1_MINUS,ǰ���ͱ����ĻҶȾ�ֵ��
%ͬʱ�����������ڸ����������ǰ�����ظ���CH1_AREA�ͱ������ظ���CH1_BKD_AREA

%CH1_MEAN_MINUS=zeros(hnum-1,vnum-1);%ǰ���ͱ����ĻҶȾ�ֵ��
CH1_RAW=zeros(hnum-1,vnum-1);%ǰ���������صĻҶ�ֵ��
CH1_BKD_RAW=zeros(hnum-1,vnum-1);%�����������صĻҶ�ֵ��
CH1_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ ǰ�����ص����
CH1_BKD_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ �������ص����
CH1_MEAN=zeros(hnum-1,vnum-1);%ǰ���ĻҶȾ�ֵ
% CH2_MEAN=zeros(hnum-1,vnum-1);
CH1_BKD_MEAN=zeros(hnum-1,vnum-1);%�����ĻҶȾ�ֵ

for i=1:hnum-1
    for j=1:vnum-1
  for ii=para2(i):para2(i+1)
    for jj=para3(j):para3(j+1)
%         for ii=para4(2*i)-t+1:para4(2*i+1)+t-1
%            for jj=para5(2*j)-t+1:para5(2*j+1)+t-1
     
            if (III(ii,jj)==1)
               CH1_RAW(i,j)=CH1_RAW(i,j)+PP1(ii,jj);%ǰ���ĻҶ�ֵ��
               CH1_AREA(i,j)=CH1_AREA(i,j)+1;%ǰ�������ظ���
            else
            
            CH1_BKD_RAW(i,j)=CH1_BKD_RAW(i,j)+PP1(ii,jj);%�����ĻҶ�ֵ��
            CH1_BKD_AREA(i,j)=CH1_BKD_AREA(i,j)+1;%���������ظ���
            end
      
          end
        end
if CH1_BKD_AREA(i,j)==0    
   CH1_BKD_MEAN(i,j)=0;%�󱳾��ĻҶȾ�ֵ
else
   CH1_BKD_MEAN(i,j)=CH1_BKD_RAW(i,j)/CH1_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ  
end
% CH1_BKD_MEAN(i,j)=CH1_BKD_RAW(i,j)/CH1_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ
if CH1_AREA(i,j)==0    
   CH1_RAW(i,j)=1; 
   CH1_MEAN(i,j)=0;%��ǰ���ĻҶȾ�ֵ
else
   CH1_MEAN(i,j)=CH1_RAW(i,j)/CH1_AREA(i,j);%��ǰ���ĻҶȾ�ֵ
  
end
CH1_MEAN(i,j)=CH1_MEAN(i,j)-CH1_BKD_MEAN(i,j);
    end
end
% xlswrite('kmeans1_CH1_AREA.xls',CH1_AREA)
% xlswrite('kmeans1_CH1_BKD_AREA.xls',CH1_BKD_AREA)
% xlswrite('kmeans1_CH1_MEAN.xls',CH1_MEAN)
% xlswrite('kmeans1_CH2_MEAN.xls',CH2_MEAN)
% xlswrite('kmeans1_CH1_BKD_MEAN.xls',CH1_BKD_MEAN)
% xlswrite('kmeans1_CH1_RAW.xls',CH1_RAW)
% xlswrite('kmeans1_CH1_BKD_RAW.xls',CH1_BKD_RAW)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ɹ���
flag=0;
medi1=median(median(CH1_AREA));
for i=1:hnum-1
    for j=1:vnum-1
        if CH1_AREA(i,j)>medi1*0.4&&CH1_AREA(i,j)<medi1*1.5
            flag=flag+1;
        end
    end
end
success1=flag/(i*j)
% success1

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_kmeans1=cputime;
%ֻ������ֵ��Ϊ������������
%�任��Сsize350*366
I02=imread(p2);
[Hh,Ww]=size(I02);

%�Աȶ���ǿ
I2=ctrenhance(p2);
%I2=double(I2);
x22=mat2gray(I2);%ת��Ϊ�Ҷ�ͼ�� 0-1
[K12,K22]=size(x22);

%%%%%%%%%%%%%%%%��ֵ�˲��õ�y1%%%%%%%%%%%%%%%%%%
m=3;
I21=medfilt2(x22,[m,m]); 
%��I2��������λ��para21����һЩ��������morphspotgrid1��ע��
%para22Ϊ����ˮƽ������λ�õľ���
%para23Ϊ���ش�ֱ������λ�õľ���
%hnum2=para21(5);%ˮƽ������
%vnum2=para21(6);%��ֱ������
%[para21,para22,para23]=morphspotgrid1(I);

[para1,para2,para3,para4,para5]=morphspotgrid(I21);
para2=round(para2);
if para2(1)==0
    para2(1)=1;
end
para3=round(para3);
if para3(1)==0
    para3(1)=1;
end
hnum=size(para2,2);%ˮƽ������
vnum=size(para3,2);%��ֱ������
% para4
% para5


I2=uint8(I21.*256);
[h,w]=size(I2);
h2=size(para4,2);
v2=size(para5,2);

t_kmeans01=cputime;
% t=3;

% hnum=hnum-3;
% h2=h2-5;
% para4=para4(1:h2);

for j=1:h
    if j<para4(2)-t
        I2(j,:)=0;
    end
    for i=3:2:h2-1
       if j>para4(i)+t&&j<para4(i+1)-t
         I2(j,:)=0;
       end
    end
    if j>para4(h2)+t
        I2(j,:)=0;
    end
end


for j=1:w
    if j<para5(2)-t
        I2(:,j)=0;
    end
    for i=1:2:v2-1
      if j>para5(i)+t&&j<para5(i+1)-t
        I2(:,j)=0;
      end
    end
    if j>para5(v2)+t
        I2(:,j)=0;
    end
end

III21=zeros(h,w);
for i1=1:hnum-1
    for j1=1:vnum-1
        for i2=1:para4(2*i1+1)-para4(2*i1)+2*t-1
            for j2=1:para5(2*j1+1)-para5(2*j1)+2*t-1
                I12(i2,j2)=I2(i2+para4(2*i1)-t,j2+para5(2*j1)-t);   
%         for i2=1:para2(i1+1)-para2(i1)+1  
%             for j2=1:para3(j1+1)-para3(j1)+1
%                 I12(i2,j2)=I(i2+para2(i1)-1,j2+para3(j1)-1); 
            end
        end
        [H,W]=size(I12);

ii=0;
num=H*W;
II2=zeros(num,1);
%����������ֵ����������II
for i=1:H
    for j=1:W
        ii=ii+1;
        II2(ii,1)=I12(i,j);
    end
end
clear I12


%ʹ��matlab�Դ���k-means����
%[idx,C,sumD,D]=kmeans(II2,2,'Replicates',10,'Display','final','MaxIter',1000,'EmptyAction', 'singleton','OnlinePhase','off');
% [idx2,C2,sumD2,D2]=kmeans(II2,2,'start','uniform');
idx2=kmeans2(II2,2,20);
% idx22=Moving_kmeans(II2,2);
idx2=2-idx2;
% idx22=2-idx22;

%����kmeans�õ����б����idx�б�ǰ���������
iii=0;
for i3=para4(2*i1)-t+1:para4(2*i1+1)+t-1
    for j3=para5(2*j1)-t+1:para5(2*j1+1)+t-1
% for i3=para2(i1):para2(i1+1)
%    for j3=para3(j1):para3(j1+1)
        iii=iii+1;
%         if idx2(1,1)==1
            if idx2(iii)==1
              III21(i3,j3)=0;
            else III21(i3,j3)=1;
            end
%         else
%             if idx2(iii)==1
%               III21(i3,j3)=1;
%             else III21(i3,j3)=0;
%             end
%         end
    end
end
    end
end

% figure,imshow(III21)
% imwrite(III21,'kmeans_cy2.tif','tif');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�󲹳�����ԭͼ�Ͻ��ָ��Ե����
%kmeans����󻮷ֱ�Ե
kmeans_edge5=edge(III21,'canny');

PP5=imread(p2);
[h1,w1]=size(PP5);
figure,imshow(PP5,[]);
for i=1:h1
  for j=1:w1
         if kmeans_edge5(i,j)==1
            hold on
         plot(j,i,'r');
         end
   end
end
% title('kmeans edge of p2 ');
time_kmean1=cputime-t_kmeans1
time_kmean01=cputime-t_kmeans01
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%�������յõ��ķָ�ͼ������������ݵĸ�������
%������������ڵĸ�������
        
PP2=imread(p2);
PP2=double(PP2);
PP2=PP2./256;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ǰ���������صĻҶ�ֵ��CH2_RAW�ͱ����������صĻҶ�ֵ��CH2_BKD_RAW
%��ǰ���������ҶȾ�ֵ����CH2_MEAN-CH2_BKD_MEAN=CH2_MINUS,ǰ���ͱ����ĻҶȾ�ֵ��
%ͬʱ�����������ڸ����������ǰ�����ظ���CH2_AREA�ͱ������ظ���CH2_BKD_AREA

%CH2_MEAN_MINUS=zeros(hnum-1,vnum-1);%ǰ���ͱ����ĻҶȾ�ֵ��
CH2_RAW=zeros(hnum-1,vnum-1);%ǰ���������صĻҶ�ֵ��
CH2_BKD_RAW=zeros(hnum-1,vnum-1);%�����������صĻҶ�ֵ��
CH2_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ ǰ�����ص����
CH2_BKD_AREA=zeros(hnum-1,vnum-1);%��¼(i,j)������ �������ص����
CH2_MEAN=zeros(hnum-1,vnum-1);%ǰ���ĻҶȾ�ֵ
CH2_BKD_MEAN=zeros(hnum-1,vnum-1);%�����ĻҶȾ�ֵ

for i=1:hnum-1
    for j=1:vnum-1
for ii=para2(i):para2(i+1)
  for jj=para3(j):para3(j+1)
%        for ii=para4(2*i)-t+1:para4(2*i+1)+t-1
%           for jj=para5(2*j)-t+1:para5(2*j+1)+t-1
              
            if (III21(ii,jj)==1)
               CH2_RAW(i,j)=CH2_RAW(i,j)+PP2(ii,jj);%ǰ���ĻҶ�ֵ��
               CH2_AREA(i,j)=CH2_AREA(i,j)+1;%ǰ�������ظ���
            else         
            CH2_BKD_RAW(i,j)=CH2_BKD_RAW(i,j)+PP2(ii,jj);%�����ĻҶ�ֵ��
            CH2_BKD_AREA(i,j)=CH2_BKD_AREA(i,j)+1;%���������ظ���
            end
            
          end
       end
if CH2_BKD_AREA(i,j)==0    
   CH2_BKD_MEAN(i,j)=0;%�󱳾��ĻҶȾ�ֵ
else
   CH2_BKD_MEAN(i,j)=CH2_BKD_RAW(i,j)/CH2_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ  
end
% CH2_BKD_MEAN(i,j)=CH1_BKD_RAW(i,j)/CH2_BKD_AREA(i,j);%�󱳾��ĻҶȾ�ֵ
if CH2_AREA(i,j)==0  
   CH2_RAW(i,j)=1; 
   CH2_MEAN(i,j)=0;%��ǰ���ĻҶȾ�ֵ
        
else
   CH2_MEAN(i,j)=CH2_RAW(i,j)/CH2_AREA(i,j);%��ǰ���ĻҶȾ�ֵ
end
CH2_MEAN(i,j)=CH2_MEAN(i,j)-CH2_BKD_MEAN(i,j);
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��ɹ���
flag=0;
medi2=median(median(CH2_AREA));
for i=1:hnum-1
    for j=1:vnum-1
        if CH2_AREA(i,j)>medi2*0.4&&CH2_AREA(i,j)<medi2*1.5
            flag=flag+1;
        end
    end
end
success2=flag/(i*j)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ʹ�þ���Ĵ�С��ͬ
% [H1,W1]=size(CH1_RAW);
% [H2,W2]=size(CH2_RAW);
% if H1>=H2
%     CH1_RAW=CH1_RAW(1:H2,1:W1);
% else  CH2_RAW=CH2_RAW(1:H1,1:W2);
% end
% [H1,W1]=size(CH1_RAW);
% [H2,W2]=size(CH2_RAW);
% if W1>=W2
%     CH1_RAW=CH1_RAW(1:H1,1:W2);
% else  CH2_RAW=CH2_RAW(1:H2,1:W1);
% end

%��log2ͼ��������ֵͼ
% kmeans_VALUE=log2((CH2_RAW)./(CH1_RAW));
kmeans_VALUE=zeros((hnum-1)*(vnum-1));
for i=1:hnum-1
    for j=1:vnum-1
%         if CH1_AREA(i,j)==0||CH2_AREA(i,j)==0||abs(log2(CH2_MEAN(i,j)/CH1_MEAN(i,j)))>10
%            kmeans_VALUE(i,j)=10;
%         else
           kmeans_VALUE(i,j)=log2(CH2_MEAN(i,j)/CH1_MEAN(i,j));
%         end
    end
end
[HH,WW]=size(kmeans_VALUE);
%�������ֵ��0�ķ����ͣ�������SVM_VAR����
%�������ֵ��0��ƫ���ͣ�������SVM_DEV����
sum1=0;
sum2=0;
for i=1:HH
    for j=1:WW
        X=[kmeans_VALUE(i,j),0];
        sum1=sum1+(std(X))^2;
        sum2=sum2+abs(kmeans_VALUE(i,j));
    end
end

kmeans_VAR=sum1
kmeans_DEV=sum2

figure;
ii=1;

for i=1:(hnum-1)
    for j=1:(vnum-1)
       hold on    
        plot(ii,kmeans_VALUE(i,j),'b*');
%         title('K-means���෨���ñ���ֵ');
        ii=ii+1;
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%��������ȡ��ı������ݶ�תΪexcel��ʽ������Ա�
xlswrite('kmeans1_CH1_AREA.xls',CH1_AREA)
xlswrite('kmeans1_CH1_BKD_AREA.xls',CH1_BKD_AREA)
xlswrite('kmeans1_CH1_MEAN.xls',CH1_MEAN)
xlswrite('kmeans1_CH1_BKD_MEAN.xls',CH1_BKD_MEAN)
xlswrite('kmeans1_CH1_RAW.xls',CH1_RAW)
xlswrite('kmeans1_CH1_BKD_RAW.xls',CH1_BKD_RAW)

xlswrite('kmeans1_CH2_AREA.xls',CH2_AREA)
xlswrite('kmeans1_CH2_BKD_AREA.xls',CH2_BKD_AREA)
xlswrite('kmeans1_CH2_MEAN.xls',CH2_MEAN)
xlswrite('kmeans1_CH2_BKD_MEAN.xls',CH2_BKD_MEAN)
xlswrite('kmeans1_CH2_RAW.xls',CH2_RAW)
xlswrite('kmeans1_CH2_BKD_RAW.xls',CH2_BKD_RAW)


xlswrite('kmeans1_VALUE.xls',kmeans_VALUE)
% xlswrite('kmeans_VAR.xls',kmeans_VAR)
% xlswrite('kmeans_DEV.xls',kmeans_DEV)
% t_kmeans1=cputime-t_kmeans1
% success1
% success2
end
